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A detailed numerical study was conducted on the dynamics and thermal response of inert, spherical panicles 
in strained, laminar, premixed hydrogetvair flames. The modeling included the solution of the steady 
conservation equations for both the gas and panicle phases along and around the stagnation streamline of an 
opposed-jet configuration, and the use of detailed descriptions of chemical kinetics and molecular transport. 
For the gas phase, the equations of mass, momentum, energy, and species are considered, while for the panicle 
phase, the model is based on conservation equations of the particle momentum balance in the axial and radial 
direction, the particle number density, and the particle thermal energy equation. The particle momentum 
equation includes the forces as induced by drag, thermophoresis, and gravity. The particle thermal energy 
equation includes the convective/conductive heat exchange between the two phases, as well as radiation 
emission and absorption by the particle. A one-point continuation method is also included in the code that 
allows for the description of turning points, typical of ignition and extinction behavior. As expected, results 
showed that the particle velocity can be substantially different than the gas phase velocity, especially in the 
presence of large temperature gradients and large strain rates. Large particles were also found to cross the gas 
stagnation plane, stagnate, and eventually reverse as a result of the opposing gas phase velocity. It was also 
shown that the particle number density varies substantially throughout the flowfield, as a result of the straining 
of the flow and the thermal expansion. Finally, for increased values of the particle number density, substantial 
flame cooling to extinction states and modification of the gas phase fluid mechanics were observed. As also 
expected, the effect of gravity was shown to be important for low convective velocities and heavy particles. 
Under such conditions, simulations indicate that the magnitude and direction of the gravitational force can 
substantially affect the profiles of the particle velocity, number density, mass flux, and temperature. © 1999 
by The Combustion Institute 


INTRODUCTION 

Combustion science has been significantly ad- 
vanced during the last 20 years mainly because 
of the evolution and extensive use of advanced 
laser diagnostics along with modeling made 
possible by supercomputers and the associated 
algorithms. Therefore, much has been learned 
in terms of the details of the underlying elemen- 
tary processes of molecular transport and chem- 
ical kinetics as well as their interaction with the 
fluid mechanics. This knowledge, however, has 
been chiefly advanced in the gas phase, in which 
the assumptions of dilute gas can be a very good 
approximation, and many gas properties can be 
accurately determined from first principles. On 
the other hand, there is a wide range of com- 
bustion phenomena of interest in which the gas 
phase interacts in many different ways with a 
liquid or solid phase. While extensive work has 
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been done on sprays, less attention has been 
given to the details of dusty reacting flows. 

Dusty flows are of particular interest for a ” 
wide range of applications. Particles can be 
present in a gas intentionally or unintentionally, 
and they can be inert or reacting. Inert particles 
in an otherwise reacting gas flow can lead to 
flame cooling and modification of the extinction 
limits of a combustible mixture. Reacting solid 
particles can release substantial amounts of heat 
upon oxidation and can be used either for 
propulsion (e.g., Al. B, and Mg) or power 
generation (coal). Furthermore, accidents can 
occur when a reacting dust accumulates in air 
and which, in the presence of an ignition source, 
can cause explosion. Such explosions can occur 
during lumber milling, in grain elevators, and in j 
mine galleries such as the one which occurred at; 
the Haswell Colliery and resulted in the famous 
Faraday and Lyell report [1]. 

Particles are also used as seeds in reacting 
flows in order to measure flow velocities by 
using either the laser Doppler velocimetry 
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< LDV) or particle image velocimetrv ( PIV) 
techniques. The particles that are used in these 
applications must be chemically inert, so that 
they do not alter the chemical composition of 
the reacting stream. Furthermore, they must be 
small and light enough to closely follow the gas 
phase, in order to assure accurate gas-phase 
„ velocity measurements. 

The addition of solid particles in a flowing gas 
stream can lead to strong couplings between the 
two phases, which can be of dynamic, thermal, 
and chemical nature. The (dynamic couplin g^? 
between the two phases is caused hy the differ, 
ent inertia possessed hv the , fluid nnd-snlid - 
phases , which results in slip between the phases 
that in turn can lead to the development of 
i nterphasia l forces and the modification of the 
velocity fields , of both phases. In addition to the^ 
phase-interaction forces, thermophoretic. cen- 
trifugal, electrostatic, magnetic, and gravita- 
tional field forces can be exerted on the parti- 
cles. Of these, electrostatic and magnetic forces 
can generally be ignored as few combustors 
possess significant electrical and magnetic fields 
and centrifugal forces can be ignored unless 
strong flow rotation exists. On the other hand, 
thermophoretic fp^ycs- which are causecLby 
STeep temperature gradients, can be important 
irl r^actlng~dusnT Hows. fiie ^ravitationa| fp^es, 
are unavoidably present and can piay~a substan- 
tial role in the overall dynamic response of the 
particle, especially for heavy panicles. Indeed, 
existing experiments on dusty flows (e.g., [2]) 
indicate that global flame properties can be 
quite different when they are conducted at 
normal- and micro -gr avity. 

To understand the dhermai effects\ >f dust on 
combustion, it must be first realized that as the 
thermaU capacity of the solid p hase is of ih e^ 
order of 1ft 3 time^.ihat of the gas phase, it will 
respond more slowly to the temperature 
changes induced by the flame and, because the 
solid and gas phases are thermally coupled, 

slows the therma l rcsponse_o£-the~gas- phase* 

Such temperature modifications can substan- 
tially affect the chemical kinetics of each phase. 
The thermal equilibration process will depend 
on the thermal properties of the two phases. 
Similarly, the emissivitv of the solid phase is 
substantially higher compared to that of the gas 


:o7 

phase so that the radiative transfer from the 
particles can result in significant energy losses. 

The ultimate challenge of this type of re- 
searcjiis to better understand the ^chemicaD 
fouphn^ between the two phases. The problem 
is complicated because the actual mechanisms 
leading to particle gasification^ c an differ de- 
pending on the chemical composition of 
particle . For exampled fcarbon burning iniriatp~ 
as a surface reaction that produces CO as a 
byproduct, which then diffuses from the surface 
and undergoes a secondary oxidation^ Metal 
powders, such as Al, B, and Mg are first "evap-~ 
orated due to elevated temperatures and their 
vapor undergoes gas phase burning with the 
surrounding [3]. The chemical coupling between 
the two phases leads to* particle size reduction. 
modificatio n. Of the gas phase species rnmpr*i. 
tjoila the elevation of the temperature of 
both phases as a result of the exothermicity of 
both the surface and gas phase reactions. How- 
ever, all particles begin this process as inert 
particles that must first be heated before they 
can themselves participate in the combustion 
process. Thus, one must first understand the 
more basic questions of the interactions of inert 
particles in a combustion environment before 
proceeding to the more complicated chemical 
effects. 

The in-detail understanding of the dynamics 
and structure of heterogeneous (sprays or par- 
ticle) flows can be only advanced by first con- 
sidering simple flow geometries, which (1) can 
be conveniently produced in the laboratory, (2) 
can be simulated with the use of detailed de- 
scription of all the all physico-chemical pro- 
cesses in both phases, and (3) are of relevance 
to practical turbulent reacting flows. Past expe- 
rience in gas phase combustion has shown that 
the opposed-jet, stagnation-type flows are 
among the most meritorious and permit an 
in-depth understanding of the details of the 
pertinent physico-chemical processes. 

A number of computational and experimen- 
tal studies on sprays and particle flows (e.g., 
[4-12]) have been conducted in stagnation-type 
configurations. Numerically, the need for a hy- 
brid fuJerian-I^grangian approach has been 
identified by Continillo and Sirignano [4], and 
the use of such approach has allowed for the 
prediction of the phenomenon of droplet flow 
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reversal [5. 6]. wh ich has been observ ed exper- 
— imentally [7. 8]. ( ^omez and Rosne x [9] have 
conducted a detailed study on the particle re- 
sponse in the opposed-jet configuration, and the 
particle thermop ho retic diff usiviugs were deter- 
mined experimentally. £ung, Law5 and co-work- 

r ers have conducted nurnencaTstudies [10. 11] 
on the effect of strain rate and temperature 


\ 
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gradients on the dynamics of inert particles, as a 
way of understanding potential errors in exper- 
imental LDV da ta that might arise fro m ther- 
mophoretih fotdes that cause the moti on of th e . 
tracei^particles’ to differ from that of the gas^ 
These studies have included the detailed de- 
scription of the particle Stokes drag and ther- 
mophoretic forces. Results have shown that 
depending on the particle size, strain rate, and 
temperature gradient, the particle velocity can 
substantially differ from the gas phase velocity, 
which indeed compromises the fidelity of exist- 
ing non-intrusive flow velocity measurement 
techniques. In these studies, the particle phase 
was not thermally coupled with the gas phase 
and was not allowed to affect the dynamics of 
the gas phase. Furthermore, the effect of gravity 
on the particle dynamics was not considered, 
since the focus of the study were the sub-micron 
sized particles, typically used as tracers, for 
which gravitational forces are negligible com- 
pared to fluid drag. 

The behavior of reacting dusty flows can also 
have a strong dynamic and thermal dependence 
on the(particle number densliyT^) which rep- 
resents TTTe number of solid particles per unit 
volume; such effects were not addressed in the 
previous numerical studies [10, 11]. Nonethe- 
less, the values of n p observed in a combustion 
environment will in most cases be very small. 
This can be easily seen by realizing that a 
stoichiometric mixture of volatile particles and 
oxidizer requires delivering masses of the same 
order of each to the flame, but since the density 
of the solid is roughly 3 orders of magnitude 
larger than that of the surrounding gas^ 
volumetric loading will be of the order of 
Consequently, it is unlikely that particles _55lL 
i nteract with one another and modeling infor- 
mation derived from studies on single particles 
in infinite fluids may be used with con fidence^. 

In view of the foregoing considerations, we 
have undertaken the task to conduct a long- 


term combined experimental and numerical 
study on the details of reacting dusty flows in 
which the effects of fluid mecnamcs. panicle 
properties, steep temperature gradients, gas 
and particle detailed kinetics, and gravity are 
systematically addressed. Given the complexity 
that reacting particles introduce, in our first 
attempt to understand such complex phenom- 
ena. the dynamic, thermal, and gravitational 
effects were independently studied for thej^ase 
of chemically inert particles; the term Gnex t? 
indicates that the particles are neither con- 
sumed nor they alter the gas phase chemistry. 
Thus in this paper, a systematic detailed numer- 
ical study is presented for inert particles embed- 
ded in strained, laminar premixed flames, which 
are stabilized in a stagnation-type flow configu- 
ration. Studies of chemically reacting particles 
will follow. 


- not beefy. 
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NUMERICAL APPROACH 


Governing Conservation Equations 


As mentioned in the Introduction, it is assumed 
that the particle number density is small enough 
so that the particles are very unlikely to encoun- 
ter one another which greatly simplifies the 
governing equations. In describing the mechan- 
ics, we will consider both the solid and the gas 
phase. The gas includes both the original oxi- 
dizer and the gaseous species generated by the 
chemical reactions. The general steady conser- 
vation equations will be presented for the case 
of chemically inert particles in axisymmetric 
coordinates. Subsequently, these equations will 
be reduced to a quasi-one-dimensional formu- 
lation appropriate for the description of stagna- 
tion-type flows. In these equations, the sub- 
scripts^^andrp’Mcorrespond to the gas and 
particleTpnase respectively. 

The overall mass continuity equation for the 
gas phase is: 


rfo/c 


i a(rp,,Vg) | *<p r »v y 
r dr dx 



The gas phase steady momentum equations 
in the axial and radial directions are. respec- 
tively: 
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( 2 ) 


(3) 


In the above, equations, u and v represent the 
velocities in the axial (x) and radial (r) direc- 
tions respectively, p. g is the mixture dynamic 
viscosity, while n £ Py a nd n./'pp are the forces 
in the axial aiuTTgcnal directions , respectively. 
exerted, upon a unit volume of gas hy. the 
particles UJj. (the expressions for the forces 
exerted oneach particle F PX and F PR in the 


axial and radial direction, respectively, will be 
given later in this section). The gravity term has 
been included only in the x-direction, which is 
parallel to the axis of symmetry of the counter- 
flow. 

The energy equation of the gas phase is given 
by 


r *Wp 


(£) 


, )T, 

tr w, V 




aTg 

dx 


d I dT *\ 

-7r[ rk >l?) 


2 VA 


dT ^ ( 37 ^ ^ 

— * + rp g 2 -jf + r 2 *k^k“>k + Qg.rad 

k«l 


In the above equ'ation. 



(4) 


is the gas phase 


mixture conductivity, c p is the mixture specific 
heat at constant pressure, c pk is the specific heat 
at constant pressure of species k , Y k is the mass 
fraction of species k, W k is the molecular weight 
of species and V'kr is the diffusion 

velocity of species k in the x and r direction 
respectively, h k is the specific enthalpy of for- 
mation of species /c, is the molar rate of 
production/destruction of species k resulting 
from all gas-phase chemical reactions, and K the 
total number of species. The term Qg^ d repre- 
sents the gas radiative loss. The term repre- 
sents, t he heat exchange (both by convection and 
conduction ) betweenlhe par tir|g* and the g as 
phases The terms (u p - u^)F PX and (v p - 
v *) pr represent the work per unit time done 
by each particle against the F PX and F PR forces. 


respectively. For a reactive dusty flow, however, 
the heat release due to chemical reactions and 
the heat transfer because of steep temperature 
gradients dominate the contributions of the 
work terms so that (u p - and (v p - 

Vg) F* P r can be considered as neg ligible. 

The gas radiation term Q g>rad is given at the 
optically thin limit by 

C?g.«d = 4 “g°t 7 1 - O- (5) 

where a is the Stefan-Boltzmann constant, T u is 
the ambient temperature, and a, is the total 
Planck’s mean absorption coefficient of the gas. 

Finally, the heat transfer term between one 
particle and the gas, Qp, is given for low (si) 
particle Reynolds number, Re p , by (e.g., [13]): 

Q p = 4ird p \ g (T g - T p ), Ch<3 (6) 
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where d p is the particle diameter. 

The conservation equation for the gas phase 
species k is 


fp e M si ^ 'Pg‘'gl 77 ) ~ ^ 


~ T- (rp e >* k r kr .) - rH'i.a* = 0 

or 

(7) 

While the presentation ofthe.gas phase equa- 
tions above is an Euleriarfone. the presentation 
of the panicle equatiops in the following will be 
in Lagrangian termsf^As first noticed by Conti- 
nillo and Sirignano [4], this can be necessary for 
the simultaneous numerical integration of the 
gas and particle phase set of equations. 

For the momentum balance equations for the 
particles, the formulation of Sung et al. [10, 11] 
was the starting point and appropriate modifi- 
cations were introduced. In the absence of 
centrifugal, electrostatic and magnetic forces, 
the complete momentum balance for a single 
particle in the axial direction, is given in a 
Lagrangian frame of reference by: 



( 8 ) 


ones studied herein [9], The effect of diffusio- 
phoresis was estimated by using the formulation 
recommended by Gomez and Rosner [9], given 
that our calculations were conducted for lean 
H 2 /air flames. As species with very different 
molecular weights are present, the diffusio- 
phoresis will be a function of the equivalence 
ratio. <t>. It was found that its maximum value is ^ 
1.4 cm/s for the <f> = 0.57 and 0.4 cnvs for the 
leaner <t> = 0.25 hydrogen/air flames, the two 
cases studied herein. Given that these values 
constitute only a minor correction to the other 
types of velocities, the diffusiophoretic effect 
was not further considered. 

In Eq. (8),^sq^epresents the Stokes drag 
force which oweTits existence to the velocity slip 
between the two phases and is given for low 
^tepby(e.g., [10, 11]): ^ 


SDX 


_ -3-ir/i.g^p(Mp ~ a g ) 


Ct&j 


This corresponds to the Type III approximation 
of Hjelmfelt and Mockros [14] with the addition 
of a slip correction factor C, which is needed in 
order to modify the Stokes’ law for the sub-p.m 
seeding particles. The Knudsen-Weber expres- 
sion for C and for all Knudsen numbers is (e.g., 
[10,11]): 



where m p is the particle mass and is given by 
ird p 

m p = 6 Ppar (9) 




and p par is the density of the particle. 

Equation (8) does not include the terms 
/Vof /hcIuJ^ representing forces due to the pressure gradient 
7 in the fluid, due to the fluid resistance to 
accelerating sphere , and due to the dra g associ- 
ated with unsteady motion. As Sung et al. [10, 
11] have argued, these terms can be neglected 
as they are proportional to the gas density, 
which is much smaller than the particle density. 
Not wkJ^Equation (8) does not also include other 
phoretic effects such as diffusiophoresis, elec- 
trophoresis, and photophoresis [9]. Any effects 
of electrophoresis was not considered as the 
particles are not charged and no electric field is 
present. The effect of photophoresis was also 
ignored as it has been shown to have negligible 
contribution for conditions analogous to the 


C = 1 + Krt[a + 0 exp (- yiKn )], (11) 

where a, 0, and y are constants derived from 
the fitting of the Knudsen-Weber formula to the 
experimental data and are equal to 1.142, 0.558, 
and 0.999, respectively [ 15]. The partic le Knud- 
sen number is defined as jKn = 2/^^ where 
i s the mean free path of the gas mo lecmgs given 
by = Opg/fpC^ (e.g., [16, 17]) where is a 
constant equal to 0.491, and c mg is the mean 
velocity of the gas molecules. Assuming that the 
molecules follow a Maxwell distribution, c mg is 
given by (e.g., [17]): 

c mg = ^8/c BOLT2 7' g /7rm g ) l/ ~. (12) 

In Eq. 12, *boltz * s Boltzmann constant 
and m g is the reduced molecule mass. For a 
multicomponent mixture m g is given by 


1 


m a 


K 

V 

k« 1 


1 


k 


( 13 ) 
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where 




m „ i, = 


? - k N 


Avogadro 


(14) 


with iV Avogadro being the Avogadro constant. 

As mentioned earlier. Eq. 10 is applicable 
only for very low Re p numbers. It was one of the 
goals of the present study, however, to investi- 
gate cases in which the particles’ Re p is larger 
compared to the studies of Sung et al. [10. 11], 
Thus, a correction was introduced to Eq. 10 
(e.g., [18]): 


f 5 Dx°— “^ (l+O. 


where Re p is defined as 


15/?Cp 687 


), 

(15) 


V Re ^ ^ ( “p — . ( 16) 

It should be noted in passing, that the relative 
velocity between the gas and particle phase is 
the appropriate characteristic velocity to use in 
the definition of Re p . Results showed that this 
correction could play an important role on the 
particle dynamic behavior as Re p increases. 

TTieterm F^x represents the (J hermog ) 
( ^Hore iipfo rc e on a spherical particle due to gas 
phase temperature gradient V7 g . In the near- 
continuum limit, Ftpx is given from Brock [19]: 


-6TT H T )g d p C,(^ + CJCn) ^ 

F TPX = — : — — ^ £ — f 

(1 + 3CJ&»)(l +2^ + 2C t ATn) ’ 

(17) 

where Tj g = p.^p g , and C m , C„ and C t constants 
with values 1.14, 1.17, and 2.18, respectively 
[20], which assure that the fitting formula covers 
satisfactorily the entire range of Kn numbers. In 
the work of Talbot, Cheng, Schefer, and Willis 

[20] , this is justified by realizing that Eq. 17 
degenerates within 3% to the collisionless limit 
(i.e., as Kn = 2 lfjjd p —* «) given by Waldmann 

[21] , The thermophoretic force E^bx k included 
because o f its signifi cance, in reacting flows , in 
which particles with diameters of the order" of 
microns are flowing against substantial temper- 
ature gradients. Then, the thermophoretic force 


has a direction — Pointin g from the high 
Ig mperatures towards the low gas temppnmr^ 
and opposes the particle motion towards th<- 
flame (e.g., 

t he term (Fgrx represents the g ravitational* 
force on each particle and is given bv: 

f GRX = - "»p g, (18) 

where # is the prevailing gravitational accelera- 
tion. In the present investigation various values 
of g were used, in order to systematically study 
the effect of gravity. Under normal gravity 
conditions.# = +981.0 cm/s~ for cases in which 
gravity opposes the particle motion (+# case), 
and# = —981.0 cm/s“ for cases in which gravity 
favors the particle motion (-# case). Baseline 
studies were also performed at zero gravity# = 
0 (0-# case). 

For a vertically oriented stagnation flow con- 
figuration in which all flame properties are 
radially uniform, the complete momentum bal- 
ance for a single particle in the radial direction 
is given in a Lagrangian frame of reference by 



where 


^sdr ~ 


-3i 77yf p (v D - v a ) 

C 


(20) 


In Eq. 20. Re p number corrections similar to the 
ones introduced in Eq. 15 are not required, 
given that around the centerline the velocities 
v p and v g are very small and their difference 
results in low Re p . 

In order to reduce th^ svstem of the gas phas e 
equali flagjrom the gjasvmmetric'I o a jjuasi-oncS 
CffinensionaP formulation similarly to Kee et al. 
[22], the quantities n p , u p , and F PX must be 
considered as functions of x only. In the ideal 
case of a strictly planar boundary layer between 


two opposed jets, this formulation describes the 
flow throughout the spatial domain between the 
two jets. In realistic experiments, however, ob- 
taining planar ge ometry throug hout is not pos- 
sible given that £edee* effect^ result from the 
unavoidable shear, which develops between the 
counterflowing jets and the ambient stagnant or 
coflowing gas. Furthermore, the experimental 
flames always have a minor curvature around 
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t he cente rline, which is result of the (pr essur e. 
{'gradient ^ which develop in the stagnation^How 
and which can affect the flatness of nozzle exit 
velocity, u exir Under such realistic conditions, 
the quasi-one dimensional formulation is a valid 
approximation at least in the immediate vicinity 
of the centerline of the system, as the symmetry 
requirement imposes that radial gradients must 
be zero asymptotically close to the centerline. 
Furthermore. Eqs. 19 and 20 can be combined 
to 


(d^\ -37T ^p(v p - r g ) 

p M dr 1 C 

and by dividing both sides by r 

or 

G P G P = ~ 3 *y p (G p - Gg). (21) 

Equation 21 was derived based on the fact that 
the radial velocities of both the gas and the 
particle must vary linearly with r, as a result of 
a Taylor's expansion asymptotically around the 
centerline (small r). The quantities G g and G p 
have been defined as 



and 



( 22 ) 


and are functions of x only, i.e., G g = G g (x) 
and G p = G p (x). 

The analysis conducted by Kee et al. [22] 
shows that if all properties are only a function of 
x then by differentiating thex-momentum equa- 
tion with respect to r and the r-momentum 
equation with respect to x, the pressure curva- 
ture J g m (1 /r)(dP/dr) emerges as an eigen- 
value of the problem. In Ref. 22, gravity was not 
included. If the gravity term pg is included in 
the ^-momentum equation, and given that this 
term is only a function ofx, at least asymptoti- 
cally around the centerline, it can be easily 
shown that J g remains an eigenvalue of the 


problem. Thus, the reduced quasi-one-dimen- 
sional gas phase mass continuin' becomes 


2p e G c 


diPgHg) 

dr 


(23) <p3^ 


The reduced quasi-one-dimensional radial mo- 
mentum equation for the gas phase becomes: 

( dG c \ , d 

W ‘dTl + pgGg ~ _ s 


dGJ 


(24) < ^3> 

It should be noted that while Eq. 23 is identical 
to the one derived for the gas phase stagnation 
flow [22], Eq. 24 is modified through the addi- 
tion of the last term in the left-hand side, which 
accounts for the gas-particle forces. 

It is important to emphasize that the presence 
of gravity does not affect the u g velocity profile 
along the centerline. This a result of the re- 
quirement that gravity acts only in the x-direc- 
tion and that radial symmetry must prevail at 
least around the centerline. Thus, for a given set 
of boundary conditions for u cxit at the two ends 
of the finite domain, G g is determined through 
the radial momentum Eq. 24 while w g is deter- 
mined from the gas phase mass continuity Eq. 

23. Gravity, however, alters the axial variation 
of the hydrodynamic pressure through the 
dP/dx term, which has to be modified at the 
exits of the two nozzles in order to maintain the 
a cxll values imposed as boundary conditions. 

This physically compensates for the presence of 
gravitational forces. Th is analysis is only valid 
either for the ideal case of infinitely large no z- 
zles resulting in strictly planar eeometrv . or at 
least around the centerline for the more rele- 
vant case of finite-size nozzles. 

The quasi-one-dimensional energy equation 
for the gas phase becomes 



dr. 


+ Pg S ^Vpk^kx -fa + S 


k« 1 

£?g,rad ~ ^ 


( 25 ) < ? 


STRUCTURE OF DUSTY REACTING FLOWS 


213 


In Eq. 25. Q^ rd d is given by Eq. 5. while the 
particle-gas heat transfer was accounted for 
larger Re p s by using a correction factor C NU : 

r- Qp = 47 rrd p A g (r g - r p )C NU . (26) 

C NU is a function of Re p that reflects the way 
that the Nusselt number increases with Re p and 
was determined as a function of In (Re p ) by 
fitting a 4th-order polynomial to the experimen- 
tal results reported by Clift, Grace, and Weber 
[23]. This correction is essential when studying 
large particles that typically operate at Re p s 
larger than one. 

The quasi-one-dimensional conservation 
equation for the gas phase species k is given by 

/dY k \ d 

<y K y p «“*( it) + ^ - o- ( 27 ) 

The formulation of the thermal energy equa- 
tion for a single particle was based on the 
assumption that the particle temperature is 
uniform and is a function of time only. This is a 
valid assumption for small particles as long as 
the solid conductivity is much larger than that of 
the gas. Thus, the particle thermal energy equa- 
tion is given in a Lagrangian frame of reference by 

<TP ^ ^p^par( ) — Qp Q p.rad ~ C?g.rad-p 

(28) 

The term Q p is given by Eq. 26. The term G p rad 
is the radiative heat loss emitted from the 
particle surface and is given by 

Cp.™i« Vptnrj- o. (29) 

where c par is the particle specific heat, Cp is the 
particle emissivity, and/t p = is the surface 
of the (spherical) panicle. 

The term Q^^-p represents the thermal 
energy, which is emitted by the gas phase and 
eventually is absorbed by each panicle. In eval- 
uating this term, each panicle at any spatial 
location was allowed to absorb heat that is 
radiated from every point in the gas phase. The 
heat absorption was determined by dividing the 
gas into layers parallel to the flame and then 
computing the radiative flux using a configura- 
tion factor for exchange between a circular disk, 
representing the gas layer, and a sphere, repre- 
senting the panicle [24]. Subsequently, Qg^-p 


was determined by summing the absorbed en- 
ergy from all gas layers. 

The mass conservation of the particles is 
described in an Eulerian frame of reference by 


2p p G 


p 


d <Pp“ P ) _ n 
dc 


where p p = njn p . For constant m p : 


2rt p G p + 


d K“p> _ n 

dr 


(30) 


Eq. 30 can be written in a Lagrangian frame of 
reference as: C - 3 (ip 

d(nJu 0 |) ** 

“t — = 0. (31) 

ds 


2n p G p + 


where(T )is the coordinate which specifies the 
Lagrangian distance traversed hv a partirlr 
since ipjection^a s the panicle may reverse di- 
rection during its flight, it can pass through the 
same Eulerian point, x y many times and its 
properties at that point can be multiple-valued. 
However, the total Lagrangian distance s , will 
steadily increase with time. Note that the veloc- 
ity in the direction of s is the absolute value of 
the Eulerian velocity |u p |. The use of s allows for 
the tracking of n r which becomes singular wher- 
ever the particles undergo path reversal. The 
following integration of Eq. 31 results in a 
formula in which the solid properties will be 
single -valued as functions of s: 


"p = 'Vinj eX P 


"t 

J 0 



dr 


= exp 

■ exp | - 

Given that 


2G P dt 


o 

' d|u 


d$ 


^ dr 


f 1 1 djMpl .. _ f* Jj£_ 

Jo * *-J o ds |« p | 

- f W 


d In |u_| = In 



and 
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exp | 

i 

then 


-In 


- “p-mj, 


pm)' 


n„ = n 


p.inj 


P - in ) 1 
U P' 


exp - 


2G p dt . 


(32) 


where n pinj and u p in) are the particle injection 
number density and velocity, i.e.. their respec- 
tive values at the nozzle exit. 

A analogous approach for the determination 
ofn p for droplets across singular reversal points 
has also been taken successfully by Gutheil and 
Sirignano [6], who studied the droplet behavior 
by dividing the computational domain in 
“sheets,” with each sheet corresponding to a 
droplet path between reversal points. Further- 
more, a rigorous mathematical explanation has 
been provided, which also proves that this sin- 
gularity is an integrable one. 

It is essential that an important clarification is 
made on the definition of the number density 
n p , in cases where there are length scales in the 
flow direction (i.e., the flame thickness), which 
are of the same size as the particles. The 
number density is technically defined as the 
number of particles counted in a volume, in the 
limit where the volume goes to zero. Presum- 
ably, that volume cannot get smaller than a 
particle diameter, but does that mean that n p is 
meaningless over such scales? The answer is 
largely one of definition of the averaging vol- 
ume and is not a problem for the present study 
where the conditions are steady at each Eule- 
rian point. Under such a condition, the extent of 
the averaging in the direction of the flow (the 
direction for which the flame thickness is a 
relevant scale) can be defined as the flow veloc- 
ity multiplied by a short period of time and can 
take that limit as the time goes to zero. Note 
that this is not a redefinition of rt p as, because 
the system is statistically stationary, it is exactly 
the same as making a volume average. Further- 
more. it should be clear that when n p is under- 
stood in this sense, there is no change in the 
meaning of any of the mass conservation equa- 
tions presented above, which is the critical issue. 

In summary, the gas phase is governed by 
Eqs. 23, 24. 25. and 27 while the panicle phase 
is governed by Eqs. 8. 21. 28. and 32. 


Description of Singular Points Through One- 
Point Continuation 

The injection of particles in flames can result in 
substantial modification of the properties of the 
gas phase: for example, the large heat capacity 
of the panicles dictates that they heat more 
slowly than the surrounding gas. If the injected 
number density n p lJnj is large enough, this cool- 
ing may lead to flame extinction. As a result, a 
turning point behavior can be expected [25] 
under cenain conditions in a flame temperature 
vs ^ p.inj diagram. Such turning-point behavior 
cannot be determined throueh varvine n as 
the extinction point is singular. This can be 
resolved by following the continuation approach 
of Nishioka et al. [26], Thus, a one-point con- 
tinuation approach was implemented by impos- 
ing a predetermined gas-phase temperature re- 
duction at one point within the flow field, so that 
the flp.inj would become part of the solution, and 
its respective injection (initial) value was re- 
moved. The internal point was chosen to be the 
location at which the temperature has maximum 
slope, following the recommendations of Nish- 
ioka et al. [26]. In the present study the one- 
point continuation approach was applied for the 
determination of n p jnj . However, this approach 
can be easily extended to study extinction 
and/or ignition as functions of the particle tem- 
perature, diameter, or any other physical prop- 
erty of the solid phase. 

Physical and Chemical Properties 

In the present investigation, the response of 
inert A1 2 0 3 particles were studied as a model 
case in premixed. atmospheric, jean H ,/air 
flames by including detailed descriptions of 
chemical kinetics and molecular transport. 

The hydrogen kinetics were taken from the 
fjRI 2.1 kinetic mechanism [27], and the mo- 
lecular properties pertinent to the calculation of 
transport coefficients were taken from the San- 
dia Transport Package code [28], 

For the calculation of the gas phase radiative 
heat transfer, the total Planck’s mean absorp- 
tion coefficient, a g , was calculated based on the 

a H,o f° r H 2 0 is given for optically thin condi- 
tions by Tien |29] and Hubbard and Tien |30] as 
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functions of temperature. Then, the total mean 
absorption coefficient, a ? is given by 

a g = a H : c/ , H : o» 

where p H , Q is the partial pressure of H : 0. More 
details on the formulation of the gas radiation 
can be found in [31, 32]. 

The physical properties for the A1 2 0 3 parti- 
cles were taken from Rosenhow, Hartnett, and 
Ganic [33]. The particle heat capacity, c par , and 
thermal conductivity, A p , were fitted as functions 
of the particle temperature by using 5th-order 
polynomials* The density of the A1 2 0 3 was 
taken as p par = 3.97 g/cm 3 and the particle 
emissivity as € p = 0.25 [33], 

Method of Solution 

The steady conservation equations for both the 
gas and particle phases were simultaneously 
solved by modifying an existing stagnation-flow 
code [34], which had been used to solve steady 
and unsteady premixed and non-premixed 
counterflowing flames. 

The gas phase equations were solved in a 
manner similar to the Sandia flame codes by 
using a damped Newton method. In terms of 
boundary conditions, the reactant, velocity, 
temperature (always 300 K), and concentrations 
at the nozzles 1 exits were specified. Further- 
more. plug ffew conditions were assumed at the 
nozzle exits More specifically, the gas phase 
velocity gradients were set equal to zero, al- 
though relaxation of such conditions is not 
expected to alter the underlying physics of the 
problem. The code was integrated into the 
CHEMKIN [35] and Transport [28] subroutine 
packages. 

The particle conservation Eqs. 8, 28, and 31 
are first-order differential equations and can be 
marched in time, (space) by using only one 
initial (boundary) condition, respectively. The 
radial momentum Eq. 21, is algebraic so that G p 
can be computed directly without any initial or 
boundary conditions. As first noticed by Conti- 
nillo and Sirignano [4], the time-marching, La- 
grangian approach for solving Eqs. 8, 28, and 31 
is preferable given that, under certain condi- 
tions, the space-marching, Eulerian approach 
cannot handle singular behaviors in space which 


mav not be s i ngular in ti me Jhis happens when 
t he particles undergo flow reversal. This is pos- 
sible when a particle’s inertia forces it to cross 
the gas stagnation plane (GSP), and penetrate 
into a counterflowing gas stream, which eventu- 
ally forces stagnation and reversal of the parti- 
cle stream. 

For the Lagrangian solution of the particle 
equations, the velocity and temperature of the 
particles were set to be equal to these of the gas 
phase at the nozzle exits (t = 0). Relaxation of 
these assumptions is, again, not expected to 
alter the physics of the problem. The Lagrang- 
ian, time marching was subsequently conducted 
with variable time steps. It. It was found that it 
was beneficial to adjust the time step Af s o that, O' 
the spatial locations resultin g from t he Lagrange — 
ian marching would coincide w ith ^- the spati a l 
locations a s deter ntinedby_th fc Fnlerian s o lu- - 
tion proce duiefqLlhe-&as phasejhus, the need 
for spatial interpolations, which inherently in- 
troduce inaccuracies, was eliminated. 

During the Lagrangian time-stepping, the 
particle axial momentum Eq. 8 was first ad- 
vanced to a new time step so that the new u p was 
determined. During the solution of Eq. 8, an 
iterative procedure was followed which would 
converge to the correct value of u p starting from 
a reasonable initial guess. This was needed 
because of the existence of the Re p number 
term in the Stokes drag force as given by Eq. 15. 

Having determined the new u p value, the parti- 
cle thermal energy Eq. 28 was then advanced for 
the determination of the value of T p at the new 
time step. Similarly to the particle axial momen- 
tum equation, an iteration scheme was also 
needed given that the fourth power of T p is 
involved in the Q p , rad term, and a closed form 
solution cannot be directly obtained. The G p 
was determined directly at all times from the 
algebraic Eq. 21. Subsequently, the particle 
number density Eq. 31 was advanced through 
the analytical expression, Eq. 32. It should be 
noted that while the behavior of Eq. 32 is 
singular at the locations of panicle flow reversal 
at which u p = 0, its integration just before and 
just after the singular point can be conveniently 
done as G p is finite throughout the particle’s 
Lagrangian path. 


ficSttoj 


/*/ 
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RESULTS AND DISCUSSION 

The interactions between the men panicles and 
a reacting gas phase can be both dynamic and 
thermal, and they are controlled by a number of 
parameters which have to be varied indepen- 
dently. These parameters are the gas phase 
injection velocity (which controls the strain 
rate) and chemical composition, and the parti- 
cle injection velocity, diameter, and number 
density. Since the present study was conducted 
for Al 2 0 ? particles only, the physical properties 
of the particles were fixed. Otherwise, these 
properties can be also independently varied. 

The numerical simulations were conducted 
for opposed-jet. atmospheric, laminar premixed 
H 2 /air flames, with a nozzle separation distance 
of L4 cm (left nozzle located at* = -*0.7 cm 
and right nozzle at* = +0.7 cm) and for <t> = 
0.57 and 0.25. The nozzle exit velocities, u CXJt , of 
the gas and particle phases were identical for 
both nozzles and varied from as low as 14 cm/s 
to as high a s 800 cm/s. 

The twin-burner assembly was considered to 
be vertical so that the gravitational forces would 
act along the direction of the system centerline. 
Given the unidirectional nature of gravity, its 
effect on the particle dynamics was considered 
for three cases. The first case was that of ( +g) 
in which the particles were injected from the 
bottom nozzle with the gravity opposing the 
particle axial motion. The second case is that of 
(-g) in which the particles were injected from 
the top nozzle with the gravity favoring the 
particle axial motion. The third case is that of 
zero gravity (0-g). The (+g) and (-g) condi- 
tions can be easily produced in the laboratory 
under normal gravity conditions. The (0-g) con- 
ditions can and will be produced in special 
NASA facilities such as a drop tower and/or a 
parabolic-trajectory aircraft. 

In all figures which follow, the particles were 
injected from the left nozzle at a spatial location 
* = -0.7 cm with a direction from left to right. 
By injecting particles from one nozzle on ly, an 
asymmetry was imposed to the system. This was 
intentional as it allows for the possibility of 
particle penetration into the opposing side of 
the GSP which leads to a variety of dynamic and 
thermal phenomena. In the previous studies of 
^ Sung et al. [10, 11] only small particles were 


studied, which nearly follow the gas phase and 
always stagnate closely to the GSP. It should be 
finally noted, that in the present calculations 
two identical premixed flames a re symmetrically 
established, one on each side of the GSP. 

In the present study, the particle diameters 
were varied fro m 0.3 to 1 00 ixm which is the range 
of interest for our scheduled experiments. Finally, 
the injection number density at the nozzle exit. 
n p mj , was varied from values as low as 10 particles/ 
cm 3 to values high enough to cool the flame and 
eventually cause local and global extinction. 


Dynamic Effects on the Particles 


The particles' dynamic response is chiefly con- 
trolled by its momentum through Eqs. 8 and 21. 
Equation 8 can be easily transformed to an 
Eulerian frame of reference, and by dividing 
both sides by Upm p and by realizing that m p is 
proportional to (d p ) 3 it can be shown that the 
spatial variation of u p is given by 

du„ A) , , + JL. 


dx Upd 


I (u, “ «p) 


(34) 


where A x and/4 2 are combinations of various 
parameters contained in Eqs. 15 and 17. A 
careful inspection of Eq. 34 reveals that as d p 
increases the contribution of the thermo- 
phoretic force is reduced. The Stokes drag 
remains important, as a larger d p implies larger 
particle inertia, which can lead to larger velocity 
difference (u g - u p ) and which partially offsets 
the effect of the diameter in the denominator. 
Equation 34 also reveals that as u p increases, 
the thermophoretic contribution is also re- 
duced. Similarly to the previous case, an arbi- 
trarily large u p can result in large slip velocities 
(u g - u p ). Thus, the effect of thermophoresis is 
expected to be important for small particles and 
low particle velocities. A comparison between the 
Stokes drag and the gravity term reveals that the 
gravitational force will be important for large 
particles and low particle velocities, as expected. 



Stokes drag and thermophoretic effects 

The effects of Stokes drag and thermophoresis 
are demonstrated in Fig. 1. which depicts the u p 
profiles of 0.3 and 5.0 p.m panicles respectively 
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Fig. 1. Gas phase and particle axial velocities profiles for a 
<b * 0.25 HVair flame with u c „ t - 100 cm/s, n p in j = 10 
panicles/cm 3 , and (a) d p - 0.3 jmn and (b) d p = 5.0 /xm with 
and without the inclusion of thermophoresis. Simulations 
also included the contributions of Stokes drag and normal 
gravity ( +g). 

for the <t> = 0.25 flame. Results indicate that for 
d p = 0.3 fxm, the particles follow closely the gas 
phase in the hydrodynamic zone and that inside 
the flame zone a substantial velocity difference 
develops between the two phases. This differ- 
ence was found to be larger for the weaker <t> = 
0.25 flame in which the extent of the thermal 
expansion and the resulting gas phase velocities 
in this expansion zone are lower compared to 
the stronger <f> = 0.57 flame. For the larger 5.0 
jim particles, it can be seen that the increased 
inertia of the particles results in a substantial 
difference between u p and w g even in the decel- 
erating hydrodynamic zone. Calculations were 
also conducted without accounting for the effect 
of the thermophoresis and the results are also 
shown in Fig. 1. It is apparent that thermo- 
phoresis is chiefly responsible for the discrep- 
ancy between u g and u p within the thermal 
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Fig. 2. Gas phase and particle axial velocities profiles for a 
<f> » 0.57 HVair flame with n p in) - 10 particles/cm 3 , for (a) 
u extt = 400 cm/s and d p = 10, 20, and 50 jim and (b) u ejol = 
S00 cm/s and d p * 10 and 20 /xm. Simulations included the 
contributions of Stokes drag, thermophoresis, and normal 
gravity (+g). 

expansion zone and that by neglecting F TPX* u p 
closely follows u g . These observations are con- 
sistent with similar previous studies (e.g., [10, 
11 ]). 

Typically, both the 0.3 and 5.0 /xm particles 
used to generate the data in Fig. 1 possess small 
panicle inenia as is desirable for tracer panicles 
used in LDV and PIV measurements, and Fig. I 
shows that these panicles can closely follow the 
gas phase and, in panicular, that they reach zero 
velocity at the GSP. By further increasing the 
panicle inenia, which may be accomplished 
either by increasing the injection velocity and/or 
the panicle diameter, it is possible that the 
particles will penetrate the GSP of the gas phase 
and stagnate at a different location, defined as 
the panicle stagnation plane (PSP). This is 
illustrated in Fig. 2 for the <f> = 0.57 flame for 
u e3U , = 400 and 800 cm/s, respectively, and for 
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various panicle sizes. Figure 2a depicts that for 
“exu ” **00 cmys, the 10 and 20 /xm panicles 
penetrate the GSP. and stagnate within the 
region of the opposing jet emerging from the 
upper burner. This leaves them with zero veloc- 
ity in a region of reversed gas flow. Conse- 
quently, the panicles reverse direction and 
again cross the GSP where they may stagnate 
again and undergo a subsequent reversal. This 
results in an oscillating motion, but one that is 
strongly damped, so that the p ailidgs ^yeiiUially, 
St agnate pn tfag GSR This oscillatory behavior 
has been previously observed for droplets [5, 6, 
8] and panicles [7] in stagnation flow configu- 
rations. The results of Fig. 2a indicate that the 
20 Jim panicles penetrate deeper into the op- 
posing jet regime than do the 10 /xm panicles 
because of their increased inenia, while the 50 
/xm panicles have a high enough inertia to reach 
all the way to the opposing nozzle exit. The 
results of Fig. 2b for the higher u cxil - 800 cm/s 
indicate that the 20 /xm panicles reach much 
deeper into the opposing jet regime compared 
to the u ejbt = 400 cm/s case, and they undergo 
the first reversal at a location very close to the 
exit of the opposing nozzle. Note that ‘kinks” 
can be observed in the 20 /xm curves for both 
injection velocities; for the a e3dt = 400 cm/s line, 
the kink can be observed on the first reversal 
point on the right hand side of the GSP, while 
for the u cxit = 800 cm/s line, the kink can b e 
seen at the second reversal point on the left side 
of the GSP. These kinks correspond to the 
Pfli nl^jwhere the oarficles pass through rhp. 
thermal expansion zones which surround the 
flames on e ither side of the GSP. Notice that, at 
those points, the gas undergoes a rapid acceler- 
ation in response to the intense heating. Thus, if 
the panicles can pass through the flames, they 
enter into a regime of reduced gas velocity 
within which they decelerate much more slowly. 
The kinks observed in Fig. 2 are a result of this 
rapid change in deceleration rate. 

Number density effects 

The variation of the particle number density, n p , 
was also studied for a wide range of conditions. 
These studies show certain peculiarities espe- 
cially around the PSP where the axial term of 
Eq. 31 becomes singular. Our method of dealing 



Fig. 3. Profiles of scaled by its injection value panicle 
number density for a <f> = 0.57 Hyair flame with u cxil = 1 14 
cm/s, V»n, = 10 particles/cnr\ and 4 P = 0.3. 5.0. 20, 50, and 
100 /xm. Simulations, included the contributions of Stokes 
drag, thermophoresis, and normal gravity (-*-£). Tempera- 
ture profile (solid line) indicates flame location. 

with the singular behavior was discussed in the 
Numerical Approach section. 

The spatial variation of n p , scaled by n p in)y is 
shown in Fig. 3 for the <f> = 0.57 flame, u CJOI = 
114 cm/s, and particle diameters 0.3, 5.0, 20, 50, 
and 100 /xm. Note that only the region to the 
left of the GSP is shown in this figure. For the 
0.3-20 /xm particles, it can be seen that n p is 
characterized by three distinct slopes. In the 
hydrodynamic region of the flow, it is constant 
or decreases slowly. It then undergoes a rapid 
reduction within the thermal expansion region 
and subsequently, this reduction becomes 
milder. The rapid reduction in the thermal 
expansion region can be understood as the gas 
stream undergoes a rapid acceleration there, 
which induces an acceleration of the particle 
stream; the reduction is then explained as Eq. 
(32) indicates that np/n p in j varies inversely to 
u p . However, the inertia of the 50 and 100 /xm 
particles is large enough to make them immune 
to all of these effects, and their velocities seem 
to decrease with a nearly-constant slope, which 
is somewhat affected by the thermal expansion. 
The results of Fig. 3 also indicate that the slope 
of n p in the hydrodynamic zone has a non- 
monotonic dependence on d p which can be 
attributed to the competition between the G p 
and d(u p )/dx terms. More specifically, it be- 
comes steeper as d p increases from 0.3 to 20 
/xm, and it becomes milder as d p further in- 
creases to 50 and 100 /xm. 
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Fig. 4. Profiles of scaled by its injection value particle 
number density for a <6 ** 0.57 H : /air flame with u eil , = 1 U 
cm/s. rt p in| =* 10 particles/cm 3 , and d p ~ 50 jtm- Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity ( +g). Arrows indicate the direction of 
the particle motion. Temperature profile (solid line) indi- 
cates flame location. 


The n p variations can be explained through 
Eq. 32 which indicates that the rt p is affected by 
both G p and d(u p )/dc. More specifically, n p 
tends to decrease as G p increases and as the 
slope d(u p )/dr becomes less negative. Physi- 
cally, a greater G p indicates more intense radial 
transport of particles, while a less negative 
d(u p )/dx slope indicates less tendency of the 
particle phase to be ‘compressed.” For very fine 
particles, e.g., d p < 1.0 jxm, the particles follow 
the gas phase very closely, and in regions in 
which the gas phase is incompressible (constant 
density) the particle phase should be also “in- 
compressible,” namely, the n p should be con- 
stant. This is reproduced by our simulations in 
Fig. 3 for the 0.3 pm particles whose n p is 
constant except in the zone, which is character- 
ized by rapid thermal expansion. For larger 
particles that do not follow the gas phase 
closely, the particle phase can have an apparent 
“compressibility” (i.e., rt p changes) even in re- 
gimes of constant gas phase density. This is 
shown in Fig. 3 for particles equal to or larger 
than 5 jim, whose n p is reduced even in the 
hydrodynamic zone of the gas phase. 

The effect of particle reversal on n p is shown 
in Fig. 4 for the <f> = 0.57 flame with = 114 
cm/s, and d p = 50 /im, conditions which result 
in several motion reversals for the particles. It is 
apparent that the particle reversal l eads to a 


Fig. 5. Gas phase and panicle axial velocities profiles for a 
<b = 0.25 HVair flame with u„ tt ~ 14 cm/s. n piBJ = 10 
particles/cm 3 , and d p = 20 4 m for ( +g). ( 0 -g), and ( -£) 
conditions. Simulations included the contributions of Stokes 
drag and thermophoresis. 


singular behavior for n r which is observed to 
increase substantially in the vicinity of the 
PSP’s. Mathematically, this is supported by Eqs. 
30, 31. and 32. At the PSP’ s, the spatial deriva- 
tive of Up becomes large and can not be bal- 
anced by the finite values of G p , which is 
independently determined from the radial mo- 
mentum equation of the particle. The results of 
Fig. 4 also indicate that as the particles reverse. 


n p rapidly drops and that continues in each 


subsequent reversal until the n p reaches negli- 
gible values. .Note that this singular behavior is 
partially a resu l t of the assumption that particles 
do not interact with one another . Obviously, the 
number density is not truly singular but is 
physicall y limited by the maximum possible 
particle packing. 




U P =dx 






Uc 


Ufr>0.d 

ClCCUrHulcftb. 




Gravitational effects 


radi'dl. 


The discussion surrounding Eq. 34 indicates 
that the effect of gravity on the particle dynam- 
ics will be substantial for low convective veloc- 
ities and heavy particles. Simulations were con- 
ducted at (+£), ( -g ), and (0-g) for the <f> = 
0.25 flame with u e]bl = 14 cm/s and d p = 5, 10, 
20, 50, and 100 jim. Selected results for the 
particle velocity u p are shown in Figs. 5 and 6 
for d p = 20 and 100 /aid, respectively. The 
results for small (i.e., £5 nm) particles are not 
shown as they generally follow the fluid motion 
and show only a small gravitational effect, only 
noticeable in the hydrodynamic zone. The devi- 
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Fig. b. Gas phase and panicle axial velocities profiles for a 
6 = 0.25 H : /air flame with u extt = 14 cm/s. n pm| = 10 
particles/cm 3 , an dd p = 100 p.m for (+g), (0-g), and (-g) 
conditions. Simulations included the contributions of Stokes 
drag and thermophoresis. Temperature profile (solid line) 
indicates flame location. 


ations disappears in the reaction zone as the 
increased drag forces, (resulting from the higher 
gas velocities that accompany the thermal ex- 
pansion) overwhelm the gravitational forces. 
The gravitational effect becomes more signifi- 
cant for the 10 and 20 ptm particles. For exam- 
ple, Fig. 5 (20 fxm particles) depicts that ( +g) 
conditions result in a negative d(u p )/dx velocity 
gradient at the nozzle exit; the magnitude of this 
change can be expected to increase with d p as 
the gravitational force increases as ( d p ) 3 while 
the drag forces increase proportional to d p . 
Thus, the (-g) condition results in a positive 
d(u p )/dr gradient at the nozzle exit which again 
increases with d p . The results for the (0-g) 
conditions indicate that u p develops initially 
under the influence of the Stokes drag, as 
expected. 

" As also expected, the gravity effect on u p was 
found to be even more dramatic for the larger 
particle diameters of 50 and 100 (Fig. 6) jim, 
respectively. For both sizes under (+g) condi- 
tions, the particles stagnate at a short distance 
(~1 mm for the present conditions) from the 
nozzle exit, before even reach the GSP, reverse, 
and eventually re-enter the nozzle. For (-g) 
and (0-g) conditions, however, the picture is 
different. Figure 6 shows that, at (0-g) the 100 
Jim particles penetrate deep into the flowfield, 
cross the GSP, stagnate inside the opposing jet 
regime, and subsequently reverse several times 



Fig. 7. For a 4> » 0.25 FWair flame with u ent = 14 cm/s, 
n p>inj * 10 particles/cm 3 , and d p = 20 (a) spatial 

variation of the particle number density scaled by its injec- 
tion value (b) spatial variation of the product of particle 
number density times particle axial velociry scaled by its 
injection value for (+g), (0-g), and (- g ) conditions. 
Simulations included the contributions of Stokes drag and 
thermophoresis. Temperature (solid line) profile indicates 
flame location. 


until they reach an equilibrium state. Figure 6 
also depicts that the large 100 particles 
under the favorable (-g) conditions, possess 
high enough inertia to shoot through the GSP 
and reach the opposing nozzle. 

A direct consequence of the gravity effect is 
the potentially large modification in the n p 
distribution, especially for large panicles. This 
can be especially seen in Fig. 7a for the 20 jun 
particles, where the substantial discrepancy of 
the n p values between the (+g), (-g), anci 
(0-g) conditions can be easily seen. The ( +g) 
conditions result in particle accumulation at the 
vicinity of the nozzle exit, while the (-g) con- 
ditions result in a substantial reduction o f n p in 
the same area. This discrepancy was found to be 
even more apparent for the 50 and 100 ^m 
particles. 

These observations can partially explain the 
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substantial differences which have been ob- 
served for the burning rates of reacting dusty 
Hows under normal- and micro-gravity condi- 
tions (e.g., [2]). The presence of gravity can 
affect the dynamic response of the particles and 
substantially modify the spatial distribution of 
the particle mass flu*. This can be seen in Fig. 
7b, which shows the spatial variation of the 
product n p *u p , (which is directly proportional 
to the particle mass flux) scaled by its injection 
value, for ( +g), (- g ), and (0-g) conditions. 
For the case of reacting particles, this is trans- 
lated to a substantial modification of the rate at 
which (solid) fuel is supplied to the reaction 
zone and thus one can expect the burning rate 
to be significantly affected by the gravitational 
orientation. 

Thermal Effects on the Particles 

The particle thermal response is controlled by 
its thermal energy equation (Eq. 28). When 
transformed to an Eulerian frame of reference, 
it can be shown that the spatial variation of the 
particle temperature T p is given by 

d7 p ^ A 3 A x 

d* ^p^p^parPpar w p^p^parPpar 

where A 3 and A a are combinations of the 
various parameters contained in Eqs. 26, 28, 
and 29. The convective-conductive terms have 
been absorbed in ihcA 3 parameter, while the 
radiation terms have been absorbed into th cA a 
parameter. Equation 35 indicates that the spa- 
tial variation of T p is inversely proportional to 
the particle velocity, density, and specific heat. 
Thus, the higher the values of u p , p pan and c pan 
the slower is the heating rate of the particle. 

Changing the particle diameter d ? has a very 
different effect on the convective/conductive 
and the radiative thermal interactions between 
the two phases. For the convective/conductive 
interaction, the particle heating rate is inversely 
proportional to d\ indicating the physical im- 
portance of the ratio between the characteristic 
length for the heat conduction (d p ) between the 
two phases and the particle volume (—d p ), 
which is a measure of the particle mass and thus 
the particle's thermal capacity. For the radiative 
interaction, the heat transfer rate is inversely 



Fig. 8. Gas phase and particle temperatures profiles for a 6 

- 0.57 H : /air flame with u e%lt = 114 cm/s, n pm) = 10 
particles/cm 3 , and d p = 0.3, 20. 50. and 100 jim. Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity ( +g). 

proportional to d p indicating the physical im- 
portance of the ratio between the particle sur- 
face (^dp), which affects the total radiative 
energy transfer, and the particle’s thermal ca- 
pacity which is again proportional to (~~d p ). 

The effect of d p on the particle heating is 
shown in Fig. 8 for the d> = 0.57 flame with u cxit 
= 114 cm/s and d p = 0.3, 20, 50, and 100 /xm. It 
is apparent that the smaller, 0.3 p,m particles 
heat up very quickly and follow the gas phase 
temperature closely. For the larger particles, a 
hysteresis exists, which becomes larger as d p 
increases. For the 50 pm particles which un- 
dergo several reversals, it is of interest to note 
that during these reversals around the GSP 
( x = 0) the particles eventually are heated up 
to temperatures which are close to those of the 
gas phase. 

The effect of particle reversal on the particle 
heating can be seen more clearly in Fig. 9 for 
the <b = 0.57 flame, with u CXJt = 400 cmys and d p 

— 20 pm. The particles are heated slowly while 
on the left side of the GSP (.t = 0). After they 
penetrate the GSP they reach a local maximum 
in temperature, and are subsequently rapidly 
cooled to their initial injection temperature just 
as they reach the first PSP. This temperature 
reduction is quite rapid as the particles undergo 
intense convective cooling by the opposing jet 
once they penetrate the GSP. Upon reversal, 
the panicles are heated, relatively quickly, back 
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Fig. 9. Gas phase and particle temperatures profiles for a 6 
= 0.57 Hyair flame with u CM = 400 cm/s. n pm) * 10 
particles/cm 3 . and d p = 20 >im. Simulations included the 
contributions of Stokes drag, thermophoresis, and normal 
gravity ( *g). Arrows indicate the direction of the panicle 
motion. 

to temperatures close to those of the gas phase, 
and as they undergo further reversals around 
the GSP, they eventually reach a maximum 
temperature that is close to the maximum gas 
temperature. Note, that during the first reversal 
path towards the flame following the first PSP 
(x — +0.22 cm), the particles are heated faster 
as compared to when the particles first enter the 
preheat zone (x ** -0.19 cm). This occurs 
because the particle velocity fields are very 
different in these two regions (as can be seen in 
Fig. 2a). Wheri the particles first enter the 
preheat zone, (0.2 cm < x < 0 cm), their 
velocities are of the order of 150 to 250 cm/s 
while in their first direction reversal (0.12 cm < 
x < 0.22 cm) their velocities range from zero to 
approximately 120 cm/s in magnitude. Conse- 
quently, the returning particles have more time 
in contact with the heated fluid and more 
rapidly rise towards the flame temperature. 

The effect of the gas and particle injection 
velocity, u CXJV on the particle heating can be 
further seen in Fig. 10 for the 6 = 0.57 flame 
with d p = 100 jxm, and u cxst = 114, 400, and 800 
cm/s. Results show that as u cxll increases, the 
two flames approach the GSP shrinking the 
region occupied by hot products, but at the 
same time, the maximum flame temperature 
increases, as is typical of these Le < 1 flames 
(e.g., [36]). However, reducing w CX]t results in 
heating the particles to higher temperatures. 



Fig. 10. Gas phase and particle temperatures profiles for a 
d> = 0.57 Hyair flame with u ejQ1 = 114, 400, and 800 cm/s. 
n pin) = 10 particles/cm 3 , and d p * 100 ^m. Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity ( +g). 

This indicates that the correspondingly larger 
residence times of the particles in the regions of 
large gas temperature are more important in 
determining the particle temperature than an 
increase of the flame temperature. Thus, for 
very high w exil (e.g., u cxil = 800 cm/s), the 
particle temperature is only slightly modified 
from its injection value despite the larger flame 
temperatures. 

The effect of gravity on the particle heating is 
shown in Fig. 11 for large particles and small 



Fig. 11. Gas phase and panicle temperatures profiles for a 
<f> = 0.25 Hyair flame with u eM = 30 cm/s, n pmj = 10 
particles/cm 3 , and d p = 100 *im for ( +g) and (0-g) condi- 
tions. Simulations included the contributions of Stokes drag 
and thermophoresis. Arrows indicate the direction of the 
panicle motion. 
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injection velocities (</> = 0.25, d p = 100 f± m. 
and u exi[ = 30 cm/s). For these heavy 

particles stagnate well before they reach the 
GSP and reverse towards the feed nozzle exit. 
During this excursion, the particles are first 
heated only to about 500 K. as they do not have 
a chance to interact with the hotter regions of 
the flame further downstream. During the re- 
turn path the particles cool more gradually than 
I the heating generating a hysteretic behavior 
that reflects the particles' large thermal capac- 
ity, which delays both their initial heating and 
the subsequent cooling. At (0-g), the particles 
penetrate the GSP and undergo several rever- 
sals before settling on the GSP. During these 
reversals, the particles follow a complex temper- 
ature variation as they sequentially interact with 
j hotter and colder gas layers. These results un- 
ambiguously show that gravity can have a po- 
tentially strong effect on the particle thermal 
x response as well as the particle dynamics. 

The relative importance of the various terms 
of Eq. 28 on the particle thermal state was also 
analyzed in detail. It was found that the convec- 
tive/conductive term, Q pi strongly dominates 
within the flame zone for all cases studied. 
Away from the flame zone, the Q p term was 
found to be comparable in magnitude to 
Qg.rad-p as radiation provides a mechanism for 
long-range interactions with the high tempera- 
ture regions of the flow. The panicle radiation 
term, (? p rad , was found to be maximum at the 
vicinity of maximum T py as it is expected, but in 
such regions its contribution is substantially 
lower than that of Q p . 

Dynamic and Thermal Effects on the Gas 
Phase 

For low panicle number density n p , the dynam- 
ics and thermal response of the gas phase are 
nearly unaffected by the presence of panicles. 
This was the case for all the results which were 
shown in Figs. 1 through 11. For these simula- 
tions, the low value of n p in ^ = 10 particles/cm 3 
was chosen in order to isolate the particle 
effects from the gas-phase effects. However, as 
n p increases, the gas phase can be affected by 
the presence of the particles both dynamically 
and thermally, as should be apparent from Eqs. 
24 and 25, respectively. 


The dynamic coupling results from the forces 
which develop between the two phases and thus, 
can affect the momentum balance of the sjas 
phase and the gas phase velocity field. This 
possibility was assessed for the <f> = 0.57 flame, 
with w exjc = 400 cm/s. d p = 50 p.m. and n p in| = 
5800 particles/cmL The simulations were con- 
ducted with and without the particle-gas force 
interaction in Eq. 28. It was found that if the 
particles penetrate the GSP. they can affect u g 
especially in a region just to the right of the 
GSP. Thus, by allowing the particle-gas interac- 
tion. u g becomes lower in magnitude. Physically, 
this is a result of the force exerted between the 
two phases, which reaches a maximum just to 
the right of the GSP where the relative velocity 
between the phases, and thus the interphasial 
drag forces, are at a maximum. As a conse- 
quence, the gas experiences a force of equal 
magnitude and opposite sign to that experi- 
enced by the particles, resulting in a reduction 
of the magnitude of w g . 

The thermal effects on the gas phase was 
found to be significant as n p increases. More 
specifically, by injecting particles at ambient 
temperature, the flames are cooled, and, as 
rtp.inj iS gradually increased, near-extinction con- 
ditions were observed. The cooling effect was 
found to be stronger for the left flame of the 
twin-flame assembly, as the left flame interacts 
first with the incoming cold particles (which are 
only injected from the left jet). It was also 
shown, that this cooling effect was different for 
particles with small and large inertia, as the 
particles with large inertia penetrate deeper 
into the flowfield and can interact thermally 
with both jets. 

The temperatures of the gas. and parti- 
cles. T p , are shown in Fig. 12 for the 6 = 0.57 
flame, with u exit = 114 cmvs. d p = 20 ^m, and 
n pm] = 10 and 24.500 particles/cm 3 . It can be 
seen that for the higher n p mjJ T g is substantially 
reduced for the left flame, while the tempera- 
ture of the right flame is nearly unaffected. For 
the higher n p inp T p is also reduced as a result of 
the reduction in flame temperature. 

In Figs. 13 and 14 results are shown for the <f> 
= 0.57 flame, with u cxit - 400 cm/s and d p = 50 
pim. For this case, the particle inertia is large 
enough that the particles can reach all the way 
to the opposing nozzle exit and thus, thermally 
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Fig. 12. Gas phase and particle temperatures profiles for a 
6 = 0.57 Hyair flame with u CM = 1 14 cm/s. n p mj ® 10 and 
24300 particles/cm 3 , and d p = 20 p.m. Simulations included 
the contributions of Stokes drag, thermophoresis, and nor- 
mal gravity (+£)• 

affect the entire flowfield. The dependence of 
the maximum flame temperature, T fmax , on 
"p.inj can be seen in Fig. 13: for all n p inj ’s, 7Y max 
is the maximum temperature of the right flame 
as the left flame undergoes more intense cool- 
ing. As expected, the results of Fig. 13 indicate 
that T (mtx is first monotonically reduced with 
n pAny ^ was found that the left flame is extin- 
guished for n p in j = 13,450 particles/cm 3 , while a 
“turning-point” behavior is observed at n p inj = 
26.627 particles/cm 3 . which indicates extinction 
of the right flame as well and as a consequence 



Fig. 13. Maximum flame temperature variation with pani- 
cle number density, for a <*> = 0.57 H 2 /air flame with u„ (l = 
400 cm/s. and d p - 50 Simulations included the 
contributions of Stokes drag, thermophorcsis. and normal 
gravity < -#). 


•900 



Fig. 14. Gas phase and particle temperature profiles for a d> 
- 0.57 FUair flame with d p = 50 ^m. u KXU = 4(X) cm/s. and 
p.mj = *0- 7,010. 13.750. and 25.258 particles/cm 3 . Simula- 
tions included the contributions of Stokes drag, thermo- 
phoresis. and normal gravity ( +£). 


global extinction. This “turning-point” extinc- 
tion behavior is typical of flames undergoing 
cooling and it could be captured by the one- 
point continuation technique that was imple- 
mented in this model. 

Figure 14 shows the thermal structures of the 
gas and particle phases for three of the cases in 
Fig. 13: for n pM = 10, 7,010, 13,760. and 25,258 
particles/cm 3 . For n p in j = 10 panicles/cm 3 case, 
T g is unaffected by the particles, while for the 
^p.inj = 7,010 particles/cm 3 case, T g is substan- 
tially reduced throughout the flame assembly 
with, as would be expected, a more profound 
effect on the left side of the flame assembly that 
first encounters the particles. T p is also smaller 
than the n p inj = 10 particles/cm 3 case, again 
reflecting the reduction in flame temperature. 
Thus, these particles remove thermal energy 
from the high-temperature regions of the gas 
phase and transfer that energy to other parts of 
the flowfield, preheating the gas flow coming 
from the right nozzle. This can be seen in Fig. 
14, as the gas phase temperature of the jet 
emerging from the right nozzle, starts increasing 
well ahead of the flame reflecting the preheat- 
ing induced by the heated particles. The n p inj = 
13,760 particles/cm 3 case describes a state just 
after the left flame has been extinguished, and 
results for both the gas and particle phases 
reflect the existence of only a single flame. The 
n pin) = 25,258 particles/cm' case describes a 
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state just before the right flame extinction, (i.e.. 
when both flames are extinguished). The sub- 
stantial reduction of the temperatures of both 
phases is apparent. 

Finally, it should be noted that for all studies 
on the effect of n p inj on the gas phase, special 
care was taken to account for the local n p 
modification when particle reversal(s) wa- 
s(were) observed. Thus, the effective n p which 
was used in the gas energy equation, was ob- 
tained by summing over the contributions to n p 
resulting by all reversal paths. 

CONCLUSIONS 

In the present study, the dynamic and thermal 
interactions between inert AJ 2 0 3 particles and 
strained, laminar, premixed H 2 /air flames were 
numerically investigated in a stagnation flow 
configuration. In this analysis, the quasi-one- 
dimensional conservation equations that de- 
scribe the ideal stagnation flow for the gas phase 
were solved. These equations include the parti- 
cle-gas phase force interaction in the momen- 
tum equation, and detailed description of chem- 
ical kinetics and molecular transport. A set of 
conservation equations was also developed for 
the particles in ideal stagnation flows. The par- 
ticle momentum balance included the effects of 
thermophoresis and gravity in addition to the 
Stokes drag force. The particle thermal energy 
balance included convective/conductive and ra- 
diative exchanges between the particle and the 
surrounding gas. A one-point continuation ap- 
proach was also implemented, allowing for the 
determination of singular, turning points, de- 
scribing flame extinction. The effects of the gas 
phase stoichiometry, strain rate, particle inertia, 
thermophoresis, and gravity were assessed. 

As expected, results show that while small, 
submicron particles follow the gas phase closely, 
larger particles fail to do so. For the small 
particles, the effect of thermophoresis was 
found to be significant in regions of large tem- 
perature gradients in the gas phase, an effect 
that was not as important for large particles. 
This is in agreement with results of previous 
studies. Particles possessing large inertia were 
found to penetrate the GSP, stagnate in the 
counterflow the opposing jet, and subsequently 


reverse direction several times around the GSP. 
in agreement with previously obser/ed flow 
reversal for droplets in stagnation-type flows. 
Finally, the results demonstrate that the param- 
eters affecting the particle velocity can also 
directly affect the particle number density which 
varies significantly throughout the flowfield. 

Simulations at ( +g) y (0-g), and ( -g) showed 
that the magnitude and direction of the gravi- 
tational force can have a strong effect on heavy 
particles moving at low velocities. In addition to 
significantly modifying the particle velocity field, 
gravity was also found to be responsible for the 
modification of the particle temperature, num- 
ber density, and mass flux distribution, which 
may partially explain experimentally observed 
effects of gravity on flame propagation in dusty 
reacting flows. 

The thermal effects between the two phases 
were also found to be strong. Results show that 
the particle heating rate is lower for high-inertia 
particles. Under conditions that result in parti- 
cle flow reversal, the particle temperature was 
shown to vary in a highly non-monotonic man- 
ner, by interacting with gas layers of vastly 
different temperatures. 

Finally, simulations with high particle number 
densities showed that the presence of the parti- 
cles can affect both the velocity and tempera- 
ture fields of the gas phase. The velocity field 
was found to be affected because of the forces 
that develop between the two phases. The effect 
on the gas phase temperature is mainly a result 
of the heat transfer between the hot flame and 
the cold particles. As the number density in- 
creases, the flames can be cooled to extinction 
states. Under conditions of high particle inertia, 
the heated particles can penetrate deep into the 
opposing jet, and for very high values of the 
number density can cause an ‘‘early’’ preheating 
of the incoming gas. 

This work way supported by NASA under grant 
NAG3-1877. Special thanks to our technical 
monitor Dr. Randall Vander WaL 
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